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ABSTRACT 

The final “giant-impact” phase of terrestrial planet formation is believed to begin 
with a large number of planetary “embryos” on nearly circular, coplanar orbits. 
Mutual gravitational interactions gradually excite their eccentricities until their 
orbits cross and they collide and merge; through this process the number of 
surviving bodies declines until the system contains a small number of planets 
on well-separated, stable orbits. In this paper we explore a simple statistical 
model for the orbit distribution of planets formed by this process, based on the 
sheared-sheet approximation and the ansatz that the planets explore uniformly 
all of the stable region of phase space. The model provides analytic predictions 
for the distribution of eccentricities and semimajor axis differences, correlations 
between orbital elements of nearby planets, and the complete N-planet distri¬ 
bution function, in terms of a single parameter, the “dynamical temperature”, 
that is determined by the planetary masses. The predicted properties are gen¬ 
erally consistent with N-body simulations of the giant-impact phase and with 
the distribution of semimajor axis differences in the Kepler catalog of extrasolar 
planets. A similar model may apply to the orbits of giant planets if these orbits 
are determined mainly by dynamical evolution after the planets have formed and 
the gas disk has disappeared. 

Subject headings: planets and satellites: dynamical evolution and stability — 
planets and satellites: formation — celestial mechanics 


1. Introduction 

It is always tempting to apply the powerful tools of statistical mechanics to macroscopic 
physical systems that exhibit some degree of regularity. The first hint of a role for statistical 
mechanics in planet-formation theory came from long-term integrations of the solar system. 
These showed that (i) the orbits of all the planets in the solar system are chaotic, with 
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Liapunov or e-folding times of a few Myr (Sussman & Wisdom 1988 Laskar 1989 Sussman 


& Wisdom 1992); (ii) the outer solar system, between Jupiter and Neptune, is “full” 


m 


the sense that there are almost no stable orbits for test particles in this region (Holman 
& Wisdom 1993; Holman 1997| ); (hi) there is a 1-2% probability that chaotic diffusion of 
Mercury’s eccentricity will lead to its loss—by ejection, collision with the Sun, or collision 
with another planet—within the next 5 Gyr (Laskar & Gastineau 2009). 


In the words of Laskar (1996), these hndings lead to the speculation that “maybe there 


was some extra planet at the early stage of formation of the solar system.. .but this led to 
so much instability that one of the planets.. .suffered a close encounter or a collision with 
the other ones. This leads eventually to the escape of the planet and the remaining system 
gets more stable. In this case, at each stage, the system should have a time of stability 
comparable with its age.” If this hypothesis is correct, the current conhguration of the solar 
system and other planetary systems might be determined, at least in part, by the statistics 
of orbital chaos. Moreover this evolution process might be approximately self-similar, in that 
the distributions of planet masses, eccentricities, inclinations, and semimajor axis differences 
in an ensemble of planetary systems would remain unchanged except for scale factors as the 
systems evolved. 

Qualitatively similar ideas have emerged in the exoplanet community, where they are 
typically called the “packed planetary systems hypothesis”. The simplest version of this 
hypothesis (Barnes & Raymond 2004 [Fang &: Margot |20i3 ) is that most planetary systems 
are “as tightly packed as possible”, that is, there is no room for additional planets on stable 
orbits. A bold prediction of this hypothesis is that if there are stable regions of phase space 
between known exoplanets then there must be undetected planets in these regions. 

Gurrent theories of terrestrial planet formation involve multiple stages and processes (for 
recent reviews see Kokubo fc Ida||2012 Morbidelli et al.||2012 Haghighipour||2013 Raymond 
et ah 2013). As the gaseous protoplanetary disk cools, dust condenses and settles into 


the disk midplane; the dust particles then accumulate into planetesimals; and the largest 
planetesimals undergo runaway growth until they dominate the gravitational scattering of 
smaller planetesimals. At this stage a phase of self-regulated or oligarchic growth begins, in 
which the largest planetesimals—now called planetary embryos—all grow at similar rates. 
The oligarchic phase ends when the reservoir of small planetesimals is exhausted. At this 
point there are typically a few dozen embryos left. The surviving embryos gradually excite 
one another’s eccentricities until their orbits cross and they collide. In this last stage— 
variously called late-stage accretion, post-oligarchic growth, or the giant-impact phase—the 
number of surviving bodies slowly declines until we are left with a small number of planets on 
well-separated, stable orbits. Laskar’s hypothesis or the packed planetary systems hypothesis 















































3 


are roughly equivalent to the assumption that the giant-impact phase tends to produce an 
ensemble of planetary systems with statistically similar properties, an idea that has been 
advanced recently from different perspectives by Malhotra (2015), Pu & Wu (2015), and 


Volk & Gladman (2015). 


The goal of this paper is to explore a simple model for the distribution of orbital elements 
in planetary systems after the giant-impact phase. The basic ansatz behind the model is that 
at the end of the giant-impact phase, the planets explore uniformly all of the stable phase 
space that is available to them. This ansatz neglects many physical processes that may play 
important roles in the late stages of planet formation (migration, gas drag, mean-motion and 
secular resonances, dynamical friction from a residual planetesimal population, hit-and-run, 
fragmenting, and catastrophic collisions, perturbations from exterior giant planets, etc.). 
The model does not describe the distribution of planetary masses following the giant-impact 
phase, only the distribution of orbits. We believe the model is useful because it yields 
clear predictions, with minimal free parameters, for several of the observables in multi¬ 
planet systems. It is, of course, no substitute for N-body simulations, but it can guide our 
interpretation of these simulations and does well at reproducing many of their results. 


The initial stages of giant planet formation are believed to be similar to those of ter¬ 
restrial planets, but the masses of giant planets are dominated by gas envelopes that are 
believed to accrete before the gaseous protoplanetary disk disappears, a few Myr after the 
birth of the star. It is possible that the distribution of giant planets continues to evolve over 
much longer times as the planets excite one another’s eccentricities. In this case the number 
of planets is whittled down mostly by ejection from the system rather than by collision as 
in the case of terrestrial planets (Chatterjee et ah 2008; Juric & Tremaine 2008 ^ Despite 


this difference, our ansatz could also apply to giant planets if this late dynamical evolution 
is the primary mechanism that determines the distribution of their orbits. 


2. The model 
2.1. The planetary system 


We shall work with a simplified model of a planetary system based on the sheared-sheet or 


Hill’s approximation 

(Spitzer & Schwarzschild||1953 

Henon |1969 , 

1970 

Petit & Henon||1986 

Binney & Tremaine 

2008 

). This simplification eliminates radial gradients in surface density. 


^The different outcomes arise because the Safronov number (eq. 
terrestrial planets. 


34) is much larger for giant planets than 
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orbital angular speed, etc., which otherwise obscure the analysis. 


In Hill’s approximation we focus on a small radial interval of the system centered on 
radius a and of width Aa. Within this interval there are N orbiting masses rrii (“planets”) 
with semimajor axes and eccentricities e*. The Hill radius of planet i is 


rm = a 



( 1 ) 


where M* is the mass of the host star. The orbital angular speed is 
We assume that the planetary orbits are coplanar or nearly so. We index the planets so 
that their semimajor axes are ordered, ai < a 2 < ■ ■ ■ < cln- The positions of the planets are 
bounded by two hxed radii Oq = a— ^Aa, otv+i = a+ |Aa, i.e., the pericenters 0^(1 —e,) > Oq 
and the apocenters aj(l + ei) < otv+i- We ignore any interactions with planets outside this 
range. 


We assume that the planet masses are small, M*, and (usually) that the number 

of planets iV 3> 1. As N grows we assume that rm Aa/N since otherwise the dynamics 
is trivial: if rui Aa/N then a system composed of planets on nearly circular orbits is 
stable, while if thi S> Aa/N there will be frequent close encounters and collisions so the 
conhguration is short-lived. This ordering condition can be satished either by letting the 
planetary masses shrink as ~ N~^ for hxed radial interval Aa (Hill’s approximation) or by 
hxing the planetary masses and letting the radial width grow as Aa ~ N (the sheared-sheet 
approximation). 


2.2. A simplified stability criterion 


The criteria for long-term stability of multi-planet systems are not completely understood. 
However, N-body integrations of systems of several planets on nearly circular and coplanar 
orbits suggest that these systems are stable if the following approximate stability criterion 


is satished ( 

Chambers et ah 

1996 

Yoshinaga et al.||1999 

Zhou et al.||2007 

Smith & Lissauer 

2009 

Funk et ah 

2010 

Pu & Wu 

2015 

): 


Oi+i - at > Acrit(t) rH-,i,i+i where rH-i,i+i = 


^i+l 


rrii + 


3M, 


mi+i\ 


1/3 




( 2 ) 


is called the mutual Hill radius. This formula assumes that a^+i — ai ai and rrii + rni+i <C 
M*. The stability parameter Acrit(t) is a dimensionless function of the age t of the system 
in units of the orbital period. Note that most of the simulations in the literature use equally 
spaced planets (either in semimajor axis or log semimajor axis) so they cannot determine 
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whether (for example) the mean semimajor axis difference or the minimnm semimajor axis 
difference determines stability. Onr stability criterion ([^ assumes that the minimnm sepa¬ 
ration shonld be used in the stability criterion, as would be natural if stability is determined 
mostly by interactions with the nearest neighbor. 


The most recent and comprehensive review of numerical determinations of A^rit is by 
Pu & Wu (2015). For typical observed exoplanetary systems, with t ~ 10^°, they estimate 


Acrit — 10.2. However, their values are somewhat below those determined by other studies 
(see Figure 3 of Pu & Wu 2015) so we shall adopt a slightly more conservative estimate 
Acrit = 11 ± 1. For N-body simulations, which typically last for t ~ 10® orbits, we use 
Acrit = 9 ± 1. 

The form of equation (|^ suggests that the minimum stable separation scales with mass 
as This empirical finding is not a consequence of rigorous dynamical arguments. For 


example, the resonance overlap criterion (Wisdom 1980; Deck et al.;2013) suggests the scaling 


however, the exponents ^ = 0.333 and | = 0.286 are sufficiently close that they are 
hard to distinguish in N-body experiments and the two scalings yield almost the same results 
for the purposes of this paper. 


For eccentric orbits, the most important combination of orbital elements that determines 
the stability of adjacent planets is the distance between the apocenter of the inner orbit and 
the pericenter of the outer one (Pu & Wu 2015 Petrovich 2015). We therefore assume that 
the system is stable if 

^ (3) 

where hi is some function of the masses ruj and ruj+i. This reduces to ([^ in the limit of 
circular orbits if 

rui + ruj+i \ 


hi Acrit (t) n 


3M* J 


( 4 ) 


An alternative criterion is given by Petrovich (2015), who has conducted extensive numerical 
experiments on the stability of systems of two planets with masses ruj/M* = 10“‘^-10“^. The 
planets are on eccentric orbits and are followed for up to 10® orbital periods. His empirical 
criterion for stability of closely spaced planets is equation ([^ with 


hi = 2.4a 


f max (rrii, mj+i) \ 






+ 0.15a. 


( 5 ) 


We normally use equation (|^ in our experiments, but we have also experimented with (|^ 
and report results with both stability criteria. 


The stability criterion (|^ implies that between any pair of adjacent planets in a stable 
system there is an excluded length hi. In a system of N planets there is a total excluded 
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length consistency with the assumptions of 12.1 we take m^v+i = mo = 0 when 

determining and hj^). Obviously if this sum exceeds Aa no stable planetary system exists. 
Thus an important parameter is the hlling factor or packing fraction 


Aa 


F = 


0<F <1. 


( 6 ) 


2.3. The distribution of orbits in phase space 

The phase space for coplanar Keplerian orbits has two degrees of freedom and actions Ji = 
I 2 = — (1 — In Hill’s approximation the actions can be taken 

to be Ji = — a) and I 2 = The canonical phase-space volume element is 

A7r^dlidl2 = de^. 

We shall assume that the planets are uniformly distributed over the available phase 
space, that is, over the phase-space volume that satisfies the stability criterion ([^. This 
assumption is reminiscent of the ergodic hypothesis, which is the basis of much of equilibrium 
statistical mechanics, but has a different interpretation in this case: it would apply, for 
example, if some process instantaneously and randomly re-distributed the planets throughout 
phase space in an ensemble of planetary systems and then only the systems having stable 
conhgurations survived. This assumption may well be incorrect: it is more likely that there 
is a slow diffusion of planets into the unstable region of phase space, which would reduce 
the phase-space density near the stability boundary. Nevertheless, the ergodic hypothesis is 
so simple and powerful that it is worthwhile to explore its implications before investigating 
more complex models with more free parameters. We shall call the models explored here 
“ergodic models”. 

Given ergodicity, the iV-planet distribution function of semimajor axes a = (ai,..., oat) 
and eccentricities e = (ei,... , cat) is 

N 

dp{a, e) = C FI{ai — dei — ao — ho) daidef — aci+i — ai — ae* — hi), (7) 

i=l 

where hf(-) is the step function, G is a normalizing constant, and we set cat+i = 0. 


2.4. Partition function 

The partition function Z is the available volume in phase space and is given by the integral 
of the right side of equation Q over the semimajor axes and eccentricities, with G = 1. 
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To evaluate this integral we first carry out the integration over semimajor axes. For given 
eccentricities the total excluded length is G = the integral over {oj} 

depends only on this total and can be written as 


cAa-G 


-Aa-G 


-Aa-G 


dxi 


dXq 


ixi 


' xn-1 


dxN = ^i7iv(Aa - G) 


( 8 ) 


where 


= 


X 


N 


if x > 0 


0 if x < 0. 

Note that Ho{x) is equal to the step function H{x). 

After doing the integral over the semimajor axes, the partition function becomes 


(9) 


Z = 


W~ 


N 




( 10 ) 


2=1 


In the sheared-sheet approximation, the integrals over van be taken from zero to infinity; 
then they can be done by induction, giving finally 


Z = 




2^{3N)\ 

where the filling factor F is defined in equation (|^. 


(1-F) 


3N 


( 11 ) 


This result can be compared to the partition function for a one-dimensional gas of N 
hard rods of length a enclosed in a box of length L (Tonks 1936 Lieb & Mattis 1966| ). In 
this case the hlling factor is F = Na/L and the partition function is 


Z 


L^(l-F)^ 

iv! 


( 12 ) 


In both cases the partition function depends only on the filling factor. The main difference 
between the gas of hard rods and our model is that the exponent N is replaced by 3N, 
reflecting the extra degrees of freedom arising from the eccentricities. 


2.5. The N-planet distribution function 

The iV-planet distribution function ([^ can be rewritten using the identity 


1 r 

H(x) = lim - / 

e-!>0+ 2TTi J 


ds 


s — le 


— exp(isa;) 


— OO 


(13) 

















(for brevity, in future equations the limit is not written explicitly but is assumed to apply 
whenever e appears). We have 


dp{si, e) = 


(3iV)!a 




dSr 


2 ( 7 ri)^+i(Aa) 3 ^(l - Sq - ie 

N 


— exp iso(ai — aei — Oq — ho) 


X 


dtti del 


dsi 


- exp [isi{ai+i - acj+i - a* - ae* - hi )], (14) 


2=1 


in which we have replaced the normalization constant C* by 1/Z so J dp{sL,e) = 1, and as 
usual e^Y+i ~ ~ ^ = n 4“ 2 An. 

The product of exponentials can be re-written as exp(i$) where 

*h 5o(ni (XCi do ho) F ^ — 1 Si(di-\-\ ne^-i-i ne^ hi) 

N 

= ^ [ai(si_i - Si) - dei{si-i -1- Si) - Sihi] + SN^d + ^Aa) - so(a - ^Aa) - soho- 
2=1 

The characteristic function of the A-planet distribution function is 


PAr(k,p)^ / dp(a,e)exp + p^ei) 


( 16 ) 


Carrying out the integrals over a*, 


^v(k,p) = — 


2iV-l(3^)!^27V poo 


7 ri^+^(Aa)^^(l — 5o — ie 

^ /*oo 

1~T I 

X 


— exp — iso{a — |Aa -|- ho) 


/ oo 

S{ki + Si_i - Si )—-^exp [-iSihi + iSiNSN{d+ |Aa)] 

-00 

roo 

X / del [*c(Pi — hSi_i — asj)]. (17) 

Jo 

Next integrate over el after multiplying by exp(—ee^) to ensure convergence: 


Pjv(k, p) = 


22^i^(3A)!a2^ dso 


27ri(Aa)^^(l — F)^^ J_^ so — ie 


exp — iso(a — ^Aa -|- ho) 


N 


X 


dsiexp[-isihi + iSiNSN{a + lAa)] 

6{ki + Si-i - Si )—- 77^73 -37- , ,32 —• (18) 


i=i 


(si - ie)(pi - asi-i - asi -f ie )2 


This result provides a complete description of the joint probability distribution of the semi¬ 
major axes and eccentricities of all N planets. 
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2.6. The one- and two-planet distribution functions 


For practical purposes, the one- or two-planet distributions are more useful than the full N- 
planet distribution. Suppose we want to study the iF-planet distribution function, starting 


with planet J + 1. Then equation (18) implies that ki = k 2 = ■ ■ ■ = kj = 0 and this implies 
that So = Si = ■ ■ ■ = sj. Moreover pi = P 2 = ■ ■ ■ = Pj = 0. Similarly kj+x+i = ■ ■ ■ = kN = 
Pjj^k+i = ■ ■ ■ = Pn = 0 and sj+k = sj+k+i = ■ ■ ■ = sn- Thus the iF-planet characteristic 
function is 


Pxikj+i, - ■ ■ ,kj+K,Pj+i, - ■ ■ Pj+k) 

2^^i^{3N)\ dsj exp [isj^Aa -a- J2j=o ^j)] 


27ri(Aa)3^(l - F)3^ 


{sj — 


X 


J+K 

n 

_i=J+l 


r°° dsiexp{-isihi)6{ki + - Sj) 

Loo (-s* - ie){pi - asi-i - asi leY 


exp [isj+K{a \ Aa - YL=. 


N 

j=J+K+l 


(19) 

h,)] 




For example let iF = 1. Then 

„ , 2i^-V3A)!a2 , , . v-J , M 

^ 7r(Aa)3^(l - F)3^ exp[^fc(a - ^Aa + Ei=oL)] 

(isexp[isAa(l — F)] 


X 


/_oo (s — — k — ie)3'^+i(p + ka — 2sa + ie^ 


( 20 ) 


with s = sj+i, k = fcj+i, p = Pj+ 1 - The 1-planet distribution function is the inverse Fourier 
transform of the characteristic function. 


pi(a,e) = 


dkdp exp[—i(/ca + pe)]Pi{k,p) 
4(3A)!a2 


(27r)2 

(3J)! [3(JV - J-^1^]! (Aa)3^'(l - “ (S " + Se)] 

X [(a + I Aa - “ eh) - a]. (21) 


If we integrate over semimajor axis, we have 


Pi(e) = / dapi{a,e) = 


12A(3A - l)a^ (1 - F - 2ea/Aa) 


3N-2 


(Ao) 


(1-F) 


3N 


F < 1. 


( 22 ) 


Notice that the result is independent of the planet number J. Moreover it is independent of 
the excluded lengths hi associated with the planet in question or its neighbors; the distribu¬ 
tion of eccentricities depends on these lengths only though the sum of the excluded lengths, 
which determines the Filing factor F through ([^. 
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Similarly, we can evaluate the two-planet distribution function for adjacent planets: 
P2(aj+i, ej+i, aj+2, ej+ 2 ) 

^_ 2^(3iV)!a^ _ , . , . 

(3J)!(3Ar-3J-6)!(Aa)3^(l-F)3^ 

X H-ij [aj+i - aej+i - (a - ^Aa -h Ei=o^i)] 

X [(a -|- I Aa — J2^=j+2^j) ~ ^J +2 ~ hej+ 2 ] 

X -ffo(oj+2 — hej+2 — oj+i — aej+i — hj+i). (23) 


We are interested in the dependence of the two-planet function on separation aj +2 —aj+i 
so we may integrate over aj+i keeping the separation fixed. We simplify the notation by 
setting aj+i -)■ a, aj +2 a', ej+i -)■ e, ej +2 e': 


p(a' — a, e, e') 


16(3A)!a^ 


7Hi{e)Hi{e')Ho \a' — a — a{e + e') — hj+i] 


(3A-5)!(Aa)3^(l-F)3^ 

X -^3^-5 [Aa(l - F) - (a' - a) - a(e + e') + /ij+i] ■ 


(24) 


Higher-order distribution functions are more complicated but could be useful; for example, 
the three-planet distribution function p(a, a', a") can be used to describe the distribution of 
semimajor axis differences a" — a if an intermediate planet is present but undetected. 


3. Comparison to simulations and observations 


We now re-cast the formulas of the preceding sections into simpler forms that can be more 
directly compared to numerical simulations, or to observations. We assume that the number 
of planets N is large, so any single system should be considered as a subsystem of one with 
a much larger radial width Aa. This approach is similar to the grand canonical ensemble 
of classical statistical mechanics. Formally we let A —)■ 00 , Aa ~ A, while the filling factor 
F ~ const. Thus, for example, a factor like (1 — x/Aa)^ —)■ exp(—Ax/Aa). 


The eccentricity distribution (eq. 22) is then 


Pi(e) = 4 exp 


Aa(l-F) 


(25) 


The mean eccentricity is (e) = 2r and we call r the dynamical temperature since it para¬ 
metrizes the level of non-circular motion in the planetary systems. The second of equations 
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(25) shows that the dynamical temperature is related to the hlling factor and the mean 


separation Aa/N; thus it is determined by the number and masses of the planets. Note also 


that the distribution (25) has a fatter tail at high eccentricities than the standard Rayleigh 


distribution, p(e) oc eexp(—ye^). 


The two-planet distribution function (24) is 


P 2 {a' — a, e, e') = Hi{e)Hi{e')HQ \a — a — a{e -\- e') — h\ 
2ar^ 

a' — a -|- a(e + e') — h 


X exp 


2ar 


(26) 


where h is the excluded length between the planets at a and a' (eq. |^. Integrating over 
semimajor axes, the joint distribution in eccentricities is 


P 2 (e,e') = ^exp 


e e' 


T 


(27) 


Thus the distribution of eccentricity of nearest neighbors is separable: the two-planet eccen¬ 


tricity distribution is the product of two one-planet distributions (eq. 25) and there is no 


correlation or anti-correlation between the eccentricities of adjacent planets. 

The distribution in semimajor axis difference and total eccentricity et = e + e' is 

poo poo 

P 2 {a' — a,et)= / ede e'de'5{et — e — e')p 2 {a'— a,e,e') 

Jo Jo 


Qar^ 


H^ia' — a — aet — h) exp 


a' — a -b aet — h 
2dr 


(28) 


Integrating over the eccentricities gives 




a' — a — h 
2dr 


where D{x) = H{x)[Qe~^ — e~‘^^{x^ + 3x‘^ -|- 6a; -|- 6)]. Taking means yields 

(a' — a) = (h) + 6ar, 

From equation ([^, for N ^ 1 the hlling factor can be estimated 

{h) 


as 


F = 


{a' — a) ’ 


(29) 


(30) 


( 31 ) 


and with this estimate equation (30) is equivalent to the second of equations (25). 
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Other authors have examined the relation between the semimajor axis differences of 
nearest neighbors and their mutual Hill radius, but typically by plotting the distribution of 
{a' — a)/h rather than a' — a — h (Fang & Margot 2013; Hansen & Murray 2013 Lissauer 


et ah 2014). In the ergodic model the second of these has a simpler interpretation: the 


distribution of a' — a — h is given by equation (29) so long as the dynamical temperature r 


is similar for all systems in the sample, whereas deriving the distribution of (a' — a)/h from 
the ergodic model requires knowledge of the distribution of h as well. 

These formulas give the distribution of eccentricities and semimajor axis differences for 
a given value of the free parameter r. Of course, different planetary systems may have 
different values of r, so fitting the formulas to a large catalog of planets is only legitimate 
if the value of r does not vary too much among the systems in the catalog. To avoid this 
difficulty, we may plot the normalized eccentricity, 

a(e' + e) 


E = 


a' — a — h’’ 


which has the distribution 


independent of r. 


p{E) = 


64^3 

OTW 


0<E <1, 


(32) 


(33) 


If the late stages of terrestrial planet formation are driven by long-term instabilities and 
collisions in all planetary systems, and all collisions result in perfect mergers (as assumed 
in most simulations), then we would expect the buildup of planets to proceed self-similarly, 
so the final filling factor E would be similar in all planetary systems. If so, then at a given 
semimajor axis, age, and stellar mass {h) ~ (eq. 0. (a' — a) ~ (eq. |^, r 
(eq. 30), and the surface density S ~ where m is the typical planet mass. 


These scalings fail when the impact velocities become sufficiently large that the collisions 
are erosive (Schlichting 2014). To explore the effect of this failure, we parametrize the 


collisional environment by the Safronov number. 


0 = 


Gm 


(34) 


where m and R are the planetary mass and radius. Here (u^) is the mean-square velocity 
relative to the local circular speed, equal to (GM*/a)(|(e^) -|- \ where (e^) and {R) are 


the mean-square eccentricity and inclination. Equation (25) gives (e^) = 6r^ so in the simple 
case where (e^) = (i^), 

e = (35) 


27 M,Rt^ 
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For the sample of Kepler planets analyzed in t3.2, {ma/{Mi,R)) ~ 0.025 so 

2 


( 0 ) = 


0.0037 




0.06 


(36) 


The typical collision becomes erosive when Vesc/^oo ^ 2-3, where Wesc = {2Gm /is the 
escape speed from the larger body and v^o = is the rms relative velocity at inhnity 

(see Leinhardt & Stewart 2012 for a mnch more thorongh analysis). Thns we expect that 


the giant-impact phase of planet formation shonld have (0)^/^ > 2-3, which in tnrn implies 
an npper limit r ;< 0.02-0.03 for the sample of Kepler planets. These argnments are ronghly 


consistent with the resnlts of Chambers (2013), who condncted N-body simnlations withont 


and with fragmentation in collisions and fonnd mean eccentricities (e) = 0.075 and 0.045 
respectively, corresponding to r = |(e) = 0.04 and 0.02. 

As discnssed in the Introdnction, the approach in this paper is to model the distribntion 
of orbital elements of a set of planets of given masses, bnt not to model how the distribntion 
of planetary masses is established. This approach is somewhat artihcial since collisions, 
merging or erosive, affect both the mass and orbit distribntions simnltaneonsly. Nevertheless, 
the conclnsion that the distribntion of planetary masses affects the orbital distribntion only 


throngh the dynamical temperatnre r or hlling factor F (the two being related by eq. 25) 
appears to be a plansible and inevitable conseqnence of the ergodic model. 


3.1. Comparison to simulations 


We compare the predictions of the ergodic model to simnlations of the late stages of terrestrial 
planet formation carried ont by Hansen & Mnrray (2013)^ They divide 2 OM 0 of material. 


distribnted between 0.05 AU and 1 AU, into ~ 30-40 planets and follow the evolntion of these 
planets for 10 Myr. Collisions are assnmed to resnlt in mergers. The calcnlation is repeated 
100 times with randomly varying initial conditions to bnild np the statistics. Following 
Hansen & Mnrray, we do not inclnde planets with a > 1.1 AU in onr analyses since these are 
mostly scattered objects that are not expected to £t the ergodic model. 


Fitting eqnation (25) to the eccentricity distribntion of 527 snrviving planets in the 


simnlations we obtain r = 0.057 ± 0.002 (l-a error, or rednction in log likelihood of ^). 


^Many other authors have also simulated the giant-impact stage of planet formation (Chambers & Wether- 


|ill|1998||Agnor et al.|1999 Chambers|2001 Kokubo et al.|200^ [Raymond et al.|2006 Morishima et al.|2008 

[Raymond et al.|2009| [Ida &: Lin|2010[ [Kokubo fc Genda|2010[ [Morishima et al.[[2010[ [Fischer fc Ciesla|2014[ 

Pfyffer et al.||2015) but the Hansen & Murray simulation offers the most direct comparison to our model. 
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Fig. 1.— The distribution of planetary eccentricities in the simulations of Hansen &: Murray 
(|2013). The curve shows the prediction of equation (25) with r = 0.057. 


With this value of r the theoretical eccentricity distribution (25) is a very good match to 


the distribution found in the simulations, as shown in Figure [T} Indeed, Hansen fc Murray 


(2013) proposed the eccentricity distribution (25) as an empirical htting formula for their 


results. 


Equation (30) relates the mean of the relative semimajor axis differences of nearest- 


neighbor pairs to r. In this case the relation depends on the stability criterion, i.e., on 
the excluded length h. Using the stability criterion ([^ with Acrit = 9 ± 1 we obtain r = 
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Fig. 2.— The distribution of semimajor axis differences (minus the excluded length) in the simula¬ 
tions of Hansen & Murray (2013). The curve shows the prediction of equation (|29l) with r = 0.067. 


0.067 ± 0.003 for the 426 pairs in the Hansen & Murray simulation; the alternative stability 
criterion (j^ yields r = 0.060. These values for the dynamical temperature are roughly 
consistent with the value determined from the eccentricities, which is an encouraging test of 
the consistency of the ergodic model. Moreover the distribution of semimajor axis differences 
is £t well by the predicted distribution (29), as shown in Figure]^ The filling factor, as 
determined by equation (31), is F = 0.30 ± 0.03 using the stability criterion (|^ and 0.36 
using 
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E 


Fig. 3.— The distribution of normalized eccentricity (eq. [3^ in the simulations of Hansen & 


Murray (2013). The curve shows the prediction of equation 


The distribution of normalized eccentricity E (eq. 32) in the simulations is shown in 
Figure]^ along with the prediction of the ergodic model. In this case the predicted distribu¬ 
tion does not match the simulations well. Part of the discrepancy is that 11% of the planet 
pairs have E > and these would be unstable according to equation]^ The other major 
discrepancy is that there are fewer planets in the simulation with E > 0.6-0.7 than the 


model predicts. These probably reflect two oversimplihcations of the ergodic model; (i) the 
stability criterion is not a sharp boundary, as assumed in the derivation; (ii) planets diffuse 
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in eccentricity towards the stability boundary, so we expect their phase-space density to be 
lower near the boundary that the uniform density predicted by the ergodic model. 



(a'-a-h)/ a 


Fig. 4.— Contours of the probability distribution for r = 0.06 of the reduced semimajor axis 
difference distribution a' — a — h and the total eccentricity (eq. 28), along with the values of these 
parameters for the nearest-neighbor pairs in the Hansen &: Murray (2013) simulation. The labels on 
each contour give the relative probability density. The region above the dashed line is not allowed 
by the stability criterion (|^ although a handful of simulated planet pairs are found there. 


Figure]^ shows contour plots of the joint probability density for the reduced semimajor 
axis difference a' — a — h and total eccentricity e -f e' of adjacent planets, derived from 
equation (28). The dynamical temperature is assumed to be r = 0.06. The red circles show 
the nearest-neighbor pairs in the Hansen & Murray (2013) simulation. In this plot, contours 
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of constant E are straight lines throngh the origin; the dashed line is -E = 1. The dehcit of 
planets at > 0.6-0.7 seen in Figure]^ is evident, since the density of red circles does not 
rise as fast as the probability density near the E = 1 line. This discrepancy is less visible 
in plots of the eccentricity or semimajor axis distribntion because these are projections onto 
the vertical or horizontal axes. 


3.2. Comparison to data 


The hnal and most important step is to compare our predictions to the actual properties of 
exoplanets. There are several obstacles to accomplishing this task: (i) Generally, eccentrici¬ 
ties are only available for planets whose orbits have been measured from radial velocities, (ii) 
The masses of most planets discovered by Kepler can only be estimated using an empirical 
mass-radius relation, and individual planet masses exhibit large deviations from the mean re¬ 
lation ( Weiss fc Marcy| 2014). (iii) There may be undiscovered planets in between the known 
members of a multi-planet system and these would affect the distribution of semimajor axis 
differences. 


Kepler planets: We have queried the NASA Exoplanet Archive for all systems containing 
more than one conhrmed planet discovered by Kepler. This sample provides 932 planets in 
362 distinct systems, containing 556 nearest-neighbor pairs. We estimate the planetary 


masses from the radii using the mass-radius relation from Weiss & Marcy (2014). We then 


compute the relative semimajor axis difference for each pair, (a' — a)/a, using a = ^(a -|- a'). 
Using the estimated masses and the stability criterion (|^ with Acrit = 11±1 we also compute 
the exclusion length h and thus determine the distribution of {a — a' — h)/a, which ought to 


be described by equation (29). 


The result is shown as the histogram in Figure The general shape of the histogram is 
similar to that of the histogram in Figure]^ from the Hansen & Murray (2013) simulations. 
One obvious difference, however, is that the histogram derived from the Kepler data has a 
signihcant number of nearest neighbor pairs with negative values of a' — a —h (48 out of 556). 
These pairs should be unstable so their presence must be explained. One possibility is that 
our stability criterion is too conservative, but this is unlikely since the analogous distribution 
from the Hansen & Murray (2013) simulations contains no pairs with a' — a — h < 0. Using 


the mass-radius relation from Fabrycky et ah (2014) or the alternative stability criterion 


(|^ does not remove this difficulty: between 31 and 63 planet pairs in the sample are still 
apparently unstable. The most plausible explanation is that the unstable pairs arise because 
of scatter in the mass-radius relation. 
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To explore this possibility we shall assume that the distribution of errors in the excluded 
length h/a (eq. is Gaussian, with standard deviation ah- Since h/a oc + for fixed 

stellar mass, an should be related to the error in roughly as au/ih) = a^i/z/ 

To estimate the latter quantity, we use the sample of 65 exoplanets with measured masses 


and radii compiled by Weiss & Marcy (2014), and compare the observed values of m}/'^ with 


the predictions from the mean mass-radius relation derived in that paper. (The mean mass 
of the Weiss & Marcy sample, {{m /= 1.68, is close to the mean for the planets in 
our sample, {{m/M^Y^'^) = 1-71, suggesting that the two samples have similar properties.) 
We hnd a^i/i /{m}/Y = 0.41. The mean excluded length in our sample is (h) = 0.25 so we 
estimate ah = 0.10 from these arguments. 


We may now convolve our predicted distribution of semimajor axis differences (29) with 


a Gaussian of width ah to find the distribution that would be obtained given realistic errors 
in the planet masses. We fit the resulting distribution to the data in Figure using the 
measurement error ah and r as fitting parameters, and find a best fit with ah = 0.12 ± 0.01, 
r = 0.028±0.004. The best-fit value for ah is remarkably close to the value obtained from the 


analysis of the Weiss & Marcy (2014) sample, ah = 0.10; thus the distribution of semimajor 


axis differences in the Kepler sample appears to be consistent with the ergodic theory once 
the scatter in the mass-radius relation is taken into account. 

The fit to the data has P^r degree of freedom is 2.7, with most of the contribution 
to xY' coming from a tail of planets with {a' — a — h)/a > 0.5 that is not present in the 
theoretical curve. The tail could arise from systems in which an intermediate planet was 
below the detection limit of the Kepler survey. 


The filling factor is F = 0.47±0.04. Using the mass-radius relation from|Fabrycky et ah 


(2014) instead of Weiss & Marcy (2014) changes this by less than the error bar, to F = 0.49. 


Using the stability criterion ([^ reduces the filling factor somewhat, to F = 0.42. 

Given the measured value r = 0.028 ± 0.004 for the Kepler semimajor axis differences. 


the mean eccentricity should be (e) = 2r = 0.05-0.06 (eq. 25), unless the eccentricities have 


been damped after the giant-impact phase is complete (Hansen & Murray 2015). So far, we 


have only limited information on the distribution of eccentricities of the Kepler planets, (i) 
Hadden & Lithwick (2014) have estimated the eccentricities using transit timing variation^ 
and obtain (e) = 0.023 ± 0.005; for planets larger than 2.5^0 the mean eccentricity is a 
factor of two smaller, or about 0.01. (ii) Moorhead et al. (2011) have estimated eccentrici¬ 
ties from the distribution of transit durations and obtain (e) = 0.1-0.25. The larger value 


^There is a typographical error in the abstract of this paper. The quantity O.OlSjlg QQ® is not the rms 
eccentricity; it is <Je which equals the rms eccentricity divided by 
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Fig. 5.— The distribution of semimajor axis differences (minus the excluded length) for confirmed 
multiple-planet systems observed by Kepler. The curve shows the prediction of equation (29) with 
r = 0.028, after convolving with a Gaussian having dispersion = 0.12 in the relative semimajor 
axis difference to account for errors in the planetary masses. 


relative to Hadden & Lithwick may arise in part because transit timing variations can only 
be measured in multi-planet systems and such systems are expected to have smaller eccen¬ 
tricities (see below); or because planets with large transit timing variations are mostly near 
strong resonances, and such planets could have a different eccentricity distribution, (iii) Sev¬ 
eral authors have estimated the mean inclination of the Kepler planets. Tremaine fc Doiig 
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(2012) find (/) < 5°; in most astrophysical disks, (e) = 1-2 times (/) so this resnlt implies 

2 



so 


find (/) ~ 1.7° and (e) ~ {/) so (e) ~ 0.03. A concern with all of these comparisons is that 
the mnlti-planet systems observed by Kepler may be biased towards low eccentricities be- 
canse eccentricity and inclination are correlated and low-inclination systems are more likely 
to have mnltiple transits. We conclnde that the observations so far are roughly consistent 
with the estimate of the mean eccentricity from the ergodic model, (e) = 2r = 0.05-0.06, 
but do not provide strong support for it. 


The hlling factor estimated from the Hansen & Murray (2013) simulations is F = 0.36± 


0.03, while the hlling factor estimated from the semimajor axis differences of Kepler planets 
is 0.4-0.5 depending on the mass-radius relation and stability criterion. If the dynamics of 
the late stages of formation of the Kepler planets are faithfully modeled by the simulations, 
we might expect the two hlling factors to be the same. The similarity of the two numbers 
is impressive but it is worthwhile to ask why they might diher. One intriguing possibility is 
that the Kepler planet masses have been overestimated; decreasing the masses by a factor 
of two would decrease the Kepler hlling factor by about 0.1 and bring it to agreement with 
the hlling factor in the simulations. 


The values of the dynamical temperature r are also similar: 0.06 in the [Hansen 


Murray (2013) simulations and 0.03 in the Kepler data. This diherence cannot be accounted 


for by diherences in the surface density or mass: the scalings described after equation (33) 
imply r ~ {m}/'^) and ((m/M 0 )^/^)HM/((^/Af©)^/^)Kepier = 1.45/1.71 = 0.85, which would 
predict that r should be 15% smaller in the simulations than in the data. A more likely 
cause of the diherence is that the simulations did not allow for fragmentation: the arguments 


following equation (36) imply that if fragmentation is present the dynamical temperature r 


cannot exceed 0.02-0.03, consistent with the Kepler data. 


Radial-velocity planets: We have queried the Exoplanets Data Explorer (^Han et ah 


2014) for multiple-planet systems discovered by radial-velocity variations in their host stars. 


This sample provides 135 planets in 55 systems or 77 nearest neighbor pairs. The distribution 
of eccentricities is ht well by equation (25) with r = 0.10; this value exceeds the upper limit 
derived after equation (36) by a factor of three or more but this is not a contradiction 
because the giant planets probably did not grow to their present masses by giant impacts. 
The distribution of semimajor axis diherences, however, is not well-ht by the ergodic model. 
There are several likely reasons for this, (i) Almost one-third of the planet pairs violate 
the stability criterion ([^ (we estimated the masses assuming the systems are edge-on, which 
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gives a lower limit); this is partly because the criterion is not valid for planets in mean-motion 
resonances, which are common among giant planets, (ii) The ergodic model is based on the 
sheared-sheet approximation, which requires that the semimajor axis difference a' — a <C a; 
this is generally not true for planets larger than a Jupiter mass since the excluded length 
h (eq. 1^ satisfies h/a = 0.96(Acrit/ll)[(?7im')/2Mjupiter]^^^- (iii) Giant planets formed 
by processes that differ from those of terrestrial planets, for example migration and the 
accretion of massive gas envelopes. As described at the end of the ergodic model is 
only expected to apply to giant planets if their orbits are mostly determined by late-stage 
dynamical evolution. 


Equation (25) suggests that in an ensemble of systems with the same hlling factor the 
mean eccentricity should vary as 1/A (more precisely, inversely with the number of planets 
per unit radius). Limbach & Turner (2015) find that for known exoplanet systems the 
eccentricity decreases with multiplicity roughly as e ~ in reasonable agreement with 

this prediction. 



Fig. 6.— As in Figure except the planet pairs plotted are Mercury-|-Venus, Venus-|-Earth, 
Earth-|-Mars. The left plot is for dynamical temperature r = 0.06 and the right for r = 0.03. 


The solar system: Figure shows contour plots of the joint probability density for the 
reduced semimajor axis difference a' — a —h and total eccentricity e -1- e' of adjacent planets, 
derived from equation (28). The dynamical temperature is assumed to be r = 0.06 in the 
left panel, the same as we found for the Kepler planets, and 0.03 on the right. The red 
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circles show the actual values for the three nearest-neighbor pairs among the four terrestrial 
planets. For r = 0.06 most of the weight of the probability density distribution lies at 
total eccentricities larger than those of the solar-system planet pairs. For r = 0.03 the 
Mercury-Venus pair has too large a semimajor axis difference relative to the probability 
density distribution. These mismatches presumably reflect the well-known difficulty that 
simulations of the giant-impact phase produce eccentricities for the terrestrial planets of the 
solar system that are too large, unless there is a residual population of small planetesimals 
to damp the eccentricities (e.g., 


Morbidelli et ah 2012). 


The outer planets of the solar system do not fit the ergodic model well, in part because 
Jupiter and Saturn are separated by only eight mutual Hill radii and so would nominally be 
unstable according to our crude stability criterion. 


Discussion 


We have described a simple model for the distribution of semimajor axes and eccentricities 
of planets. The model assumes that the last phase of terrestrial planet formation was the 
giant-impact phase, and is based on the simple ansatz that planets are uniformly distributed 
over the volume of phase space in which their orbits are stable for the lifetime of the plan¬ 
etary system (the “ergodic model”). Our model yields predictions, in terms of a single free 
parameter r that we call the dynamical temperature, for the distribution of eccentricities 


(eq. 25), the distribution of separations of nearest neighbors (eq. 29), and more generally 
for any property derivable from the complete V-planet distribution function. For example, 
the ergodic model predicts that in a given system the eccentricities should be independent 
of planetary mass. N-body simulations generally report only a weak negative correlation 
between eccentricity and mass (e.g.. 


Chambers 2013). 


The ergodic model has many limitations, (i) It does not account for the likely influence of 
giant planets at larger radii on the formation of terrestrial planets, (ii) Different planetary 
systems may have different dynamical temperatures r and fitting the data from a large 
ensemble of systems to the ergodic model will only work well if most systems have similar 
temperatures, (hi) Our analysis is based on the sheared-sheet or Hill’s approximation and 
hence does not work well when the planetary masses are large enough that the relative 
separations {a' — a)/a or the fractional excluded lengths h (eq. are of order unity or larger; 
typically this occurs for planets more massive than Jupiter, (iv) The stability criteria (|^ and 
(|^ are only approximate and probably no simple criterion perfectly separates stable from 
unstable orbits in a multi-planet system ( Petrovich|2015[ ) . (v) Our assumption that the orbits 
are distributed uniformly over the stable part of phase space is shaky, essentially because 
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the neighboring nnstable regions represent an absorbing barrier rather than a reflecting one. 
A better approximation would require following the diffusion of planet orbits through phase 
space, but implementing this would require a reliable model for the rate of this diffusion and 
how it depends on the orbits of the planet and its neighbors. 

In this paper we examine only the distribution of orbital elements and not the distribu¬ 
tion of planetary masses that emerges in the giant-impact phase, and any complete statistical 
model of this phase should predict both. An interesting question is whether the system de¬ 
scribed in §2.1[ based on the sheared-sheet or Hill’s approximation, exhibits long-range order 
in the masses as —>■ oo; in other words if the radial width Aa ~ N but the average surface 
density is fixed, does the mass of the most massive planet grow as rumax ~ and if so what 
is the critical exponent fc? 


Although the ergodic model is simple and physically plausible, there are other ap¬ 
proaches to the statistical mechanics of planet formation. An interesting alternative is due 
to Laskar (2000), who pointed out that in the secular approximation there is an important 


conserved quantity: the angular-momentum deficit. 


N 


C = - (1 - (37) 

i=l 

that is, the difference between the total angular momentum of the planets and the angular 
momentum that they would have on circular, coplanar orbits with the same semimajor 
axes. He hypothesized that the eccentricities and inclinations of planets evolve randomly, 
subject to conservation of the angular-momentum deficit, until there is a binary collision; 
in each collision the angular-momentum deficit is reduced; and collisions cease and the 
system becomes permanently stable once the angular-momentum deficit becomes too small 
to allow any more close encounters. Laskar’s model is not unique, since it requires an ad 
hoc prescription for the “random” evolution of the orbits; however, once this prescription 
is implemented it is straightforward to predict both the masses and the orbital elements of 
the planets in an ensemble of systems. Laskar’s model differs from ours in that it generally 
predicts a strong anti-correlation between eccentricity and mass. 


One-dimensional models are powerful tools in statistical mechanics because they are 
often much simpler to solve than their three-dimensional analogs (Lieb & Mattis 1966). It 
is remarkable that the giant-impact phase of planetary formation is most naturally modeled 
in one physical dimension, radius, basically because (i) the systems are nearly flat; (ii) 
orbital and apsidal motion effectively averages the orbital and collisional dynamics over the 
azimuthal angle. Although the ergodic model contains many simplifications, more general 
and accurate analyses of the statistical mechanics of the giant-impact phase could still be 
one-dimensional in this sense and therefore amenable to analytic treatments. 
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Exact models in statistical mechanics are mainly useful because they provide insight 
into the behavior of real systems, rather than because they are accurate representations of 
them. Similarly, the ergodic model presented here is mainly useful because it provides simple 
predictions for many properties of the 1-, 2-, and even iV-planet joint distribution of orbital 
elements. These predictions encapsulate some, though certainly not all, of the physics of 
the late stages of planet formation, and therefore should help to organize and interpret both 
observations of planet orbits and numerical simulations of planet formation. 

This research was initially stimulated by conversations with Renu Malhotra. I thank 
Brad Hansen for comments on the manuscript, and for providing the results from his simula¬ 
tions. I am particularly grateful to Cristobal Petrovich for discussions, insight, and pointers 
to the literature, and for carrying out exploratory N-body integrations on my behalf. 


REFERENCES 

Agnor, C. B., Canup, R. M., & Levison, H. F. 1999, Icarus, 142, 219 
Barnes, R., & Raymond, S. N. 2004, ApJ, 617, 569 

Binney, J., & Tremaine, S. 2008, Galactic Dynamics (Princeton: Princeton University Press) 

Chambers, J. E. 2001, Icarus, 152, 205 

Chambers, J. E. 2013, Icarus, 224, 43 

Chambers, J. E., & Wetherill, G. W. 1998, Icarus, 136, 304 

Chambers, J. E., Wetherill, G. W., & Boss, A. P. 1996, Icarus, 119, 261 

Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580 

Deck, K. M., Payne, M., & Holman, M. J. 2013, ApJ, 774, id. 129 

Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, id. 146 

Fang, J., & Margot, J.-L. 2012, ApJ, 761, id. 92 

Fang, J., & Margot, J.-L. 2013, ApJ, 767, id. 115 

Figueira, P., Marmier, M., Bone, G., et al. 2012, A&A, 541, id. A139 

Fischer, R. A., & Ciesla, F. J. 2014, Earth and Planetary Science Letters, 392, 28 



Funk, B., Wuchterl, G., Schwarz, R., Pilat-Lohinger, E., & Eggl, S. 2010, A&A, 516, id. A82 
Hadden, S., & Lithwick, Y. 2014, ApJ, 787, id. 80 

Haghighipour, N. 2013, Annual Review of Earth and Planetary Sciences, 41, 469 

Han, E., Wang, S. X., Wright, J. T., et al. 2014, PASP, 126, 827 

Hansen, B.M.S., & Murray, N. 2013, ApJ, 775, id. 53 

Hansen, B.M.S., & Murray, N. 2015, MNRAS, 448, 1044 

Henon, M. 1969, A&A, 1, 223 

Henon, M. 1970, A&A, 9, 24 

Holman, M. J. 1997, Nature, 387, 785 

Holman, M. J., & Wisdom, J. 1993, AJ, 105, 1987 

Ida, S., & Lin, D. N. C. 2010, ApJ, 719, 810 

Johansen, A., Davies, M. B., Church, R. P., & Holmelin, V. 2012, ApJ, 758, id. 39 
Juric, M., & Tremaine, S. 2008, ApJ, 686, 603 
Kokubo, E., & Genda, H. 2010, ApJ, 714, L21 

Kokubo, E., & Ida, S. 2012, Progress of Theoretical and Experimental Physics, 2012, id. 
01A308 

Kokubo, E., Kominami, J., & Ida, S. 2006, ApJ, 642, 1131 
Laskar, J. 1989, Nature, 338, 237 

Laskar, J. 1996, in Dynamics, Ephemerides, and Astrometry of the Solar System, lAU 
Symposium 172, ed. S. Ferraz-Mello, B. Morando, and J.-E. Arlot (Dordrecht, NL: 
Kluwer), 75 

Laskar, J. 2000, Physical Review Letters, 84, 3240 
Laskar, J., & Gastineau, M. 2009, Nature, 459, 817 
Leinhardt, Z. M., & Stewart, S. T. 2012, ApJ, 745, id. 79 



- 27- 

Lieb, E. H., & Mattis, D. C. 1966, Mathematical Physics in One Dimension (New York: 
Academic Press) 

Limbach, M. A., & Turner, E. L. 2015, PNAS, 112, 20 

Lissauer, J. J., Dawson, R. h, & Tremaine, S. 2014, Nature, 513, 336 

Malhotra, R. 2015, arXiv: 1502.05011 

Moorhead, A. V., Ford, E. B., Morehead, R. C., et ah 2011, ApJS, 197, id. 1 

Morbidelli, A., Lunine, J. I., O’Brien, D. P., Raymond, S. N., & Walsh, K. J. 2012, Annual 
Review of Earth and Planetary Sciences, 40, 251 

Morishima, R., Schmidt, M. W., Stadel, J., & Moore, B. 2008, ApJ, 685, 1247 

Morishima, R., Stadel, J., & Moore, B. 2010, Icarus, 207, 517 

Petit, J.-M., & Henon, M. 1986, Icarus, 66, 536 

Petrovich, C. 2015, in preparation 

Pfyffer, S., Alibert, Y., Benz, W., & Swoboda, D. 2015, arXiv: 1502.04260 
Pu, B., & Wu, Y. 2015, arXiv:1502.05449 

Raymond, S. N., Quinn, T., & Lunine, J. I. 2006, Icarus, 183, 265 

Raymond, S. N., O’Brien, D. P., Morbidelli, A., & Kaib, N. A. 2009, Icarus, 203, 644 

Raymond, S. N., Kokubo, E., Morbidelli, A., Morishima, R., & Walsh, K. J. 2014, in Proto¬ 

stars and Planets VI, eds. H. Beuther, R. Klessen, C. DuIIemond, T. Henning (Tucson: 
University of Arizona Press). Also arXiv: 1312.1689, 

Schlichting, H. E. 2014, ApJ, 795, id. L15 

Smith, A. W., Sz Lissauer, J. J. 2009, Icarus, 201, 381 

Spitzer, L., Jr., & Schwarzschild, M. 1953, ApJ, 118, 106 

Sussman, G. J., & Wisdom, J. 1988, Science, 241, 433 

Sussman, G. J., &: Wisdom, J. 1992, Science, 257, 56 

Tonks, L. 1966, Phys. Rev., 50, 955 





Tremaine, S., & Dong, S. 2012, AJ, 143, id. 94 
Volk, K., & Gladman, B. 2015, arXiv: 1502.06558 
Weiss, L. M., & Marcy, G. W. 2014, ApJ, 783, id. L6 
Wisdom, J. 1980, AJ, 85, 1122 

Yoshinaga, K., Kokubo, E., & Makino, J. 1999, Icarus, 139, 328 
Zhou, J.-L., Lin, D. N. G., & Sun, Y.-S. 2007, ApJ, 666, 423 


This preprint was prepared with the AAS lAILX macros v5.2. 




